Podstawy bioinformatyki w Mathematice

Ewelina Kubicz
Biofizyka III rok; studia pierwszego stopnia
Wydział Fizyki, Astronomii i Informatyki Stosowanej
Uniwersytet Jagielloński
Kraków 2011
Mathematica 8.0

Spis treści

1. Spis treści           			                                                    str. 2
2. Czym zajmuje się bioinformatyka                                                             str. 3
3. Teoria                         				                                     str. 4
4. Transkrypcja                                              	                                     str. 5
5. Translacja             				                                     str. 6
6. Przykład transkrypcji i translacji                                                                str. 7
7. Replikacja, delecja, mutacja i insercja DNA                                          str. 8 
8. Przykład 					                                     str. 9  
9. Podobieństwa sekwencji DNA                                                                 str. 10 
10. Optymalne dopasowanie dwóch sekwencji DNA                              str. 11 
11. Podsumownie 				                                    str. 12 
12. Bibliografia 					                                    str. 13

Czym zajmuje się bioinformatyka?

Bioinformatyka – jest to zastosowanie narzędzi informatycznych i matematycznych w rozwiązywaniu zagadnień z nauk biologicznych, chemicznych. Podstawowymi dziedzinami tej nauki są: 
– Genomika – zajmuje się poznaniem sekwencji i funkcjii genomu oraz zależności i reakcji zachodzących wewnątrz niego
– Proteomika – badanie struktury białek oraz jej wpływu na pełnione prze nie funkcje
– Projektowanie leków

Bioinformatyka zajmuje się również katalogowaniem i przetwarzaniem informacjii biologicznych zawartych w bazach danych takich jak:
– analiza sekwencji DNA, genomów i białek
– analiza relacji ewolucyjnych pomiędzi zbiorami sekwencji (filogenetyka)
– analiza ekspresji genów
– przewidywanie i katalogowanie struktury i funkcji genów i białek

Bioinfromatyczne bazy danych zawierają ogromne ilości różnego rodzaju informacji z dziedziny biologii, chemii. Dla przykładu:
,,Zapis genomu człowieka, to stos pudełek z CD o wysokości 750 m (ok. 50 TB)" – Świat Nauki, Październik 2000.
Główne bazy danych to:
– PDB – Protein Data Bank
– GenBank
– HGMD – Human Gene Mutation Database
 Mathematica również posiada ogromne bazy danych np.: ChemicalData, ProteinData, GenomeData.

Teoria

Aminokwasy

Za pomocą bazy ChemicalData możemy zdobyć ogromną ilość danych na temat różnych związków chemicznych takich jak:

ChemicalData["Properties"]

Przedstawiona poniżej tabela pokazuje podstawowe informacje o aminokwasach. Z uwagi, iż wszystkie dane w Mathematice są w języku angielskim dodałam również polskie nazwy.

amino= ChemicalData["AminoAcids"]

Length[amino]

nazwypolskie={"Glicyna", "(L)Alanina", "(L)Seryna", "(L)Prolina", "(L)Walina","(L)Treonina","(L)Cysteina", "(L)Izoleucyna","(L)Leucyna","(L)Asparagina","(L)Kwas asparaginowy","(L)Glutamina", "(L)Lizyna", "(L)Kwas glutaminowy","(L)Metionina", "(L)Histydyna", "(L)Fenyloalanina", "(L)Arginina", "(L)Tyrozyna", "(L)Tryptofan"};

Aminokwasy = Join[{{"Nazwa aminokwasu", "Wzór strukturalny", "Struktura przestrzenna", "Nazwa polska"}}, Table[
{amino[[i]],
ChemicalData[amino[[i]], "ColorStructureDiagram"],ChemicalData[amino[[i]], "MoleculePlot"], nazwypolskie[[i]]
} , {i,1,20,1}]]// TableForm

Kwasy nukleinowe

Są to związki organiczne zbudowane z reszty cukrowej (pentozy), reszty fosforanowej i zasady azotowej purynowej (adenina,guanina) lub pirimidynowej (cytozyna, tymina, uracyl). Podstawowe typy naturalnych kwasów nukleinowych to:
– Kwas deoksrybonukleinowy – DNA
– Kwas rybonukleinowy – RNA
Ich główną funkcją jest przechowywanie inforamcji genetycznej.

Transkrypcja

Pierwszym z procesów realizacji informacji genetycznej jest transkrypcja. Polega ona na przepisaniu informacji genetycznej z DNA na RNA. Przeprowadzenie tej reakcji umożliwia enzym polimeraza RNA.
Funkcja dokonująca transkrypcji:

transkrypcja[dna_?StringQ,n_]:=Module[{dna1,dl,dl2,rna},
dna1=StringReplace[dna,RegularExpression["[^ACGTX]"]->"X"];
dl= StringLength[dna1]-n+1;
If[dl>2,dl2=dl-Mod[dl,3];
rna=StringReplace[dna1,"T"->"U"];
Table[StringTake[rna,{i,i+2}],{i,n,n+dl2-1-2,3}],
"ERROR!"]]

transkrypcja1[dna_?StringQ]:=transkrypcja[dna,1]

Następnie sekwencje zapisujemy w postaci tzw. ,,trójek nukleotydów" czyli kodonów, przy czym stop oznacza tzw. ,,kodon stop":

DNA=Tuples[{"A", "C", "G", "T"},3];
RNA=DNA/."T"->"U";

stop="!STOP!";

codons=Table[StringJoin[RNA[[i]]],{i,1,64}]

Translacja

Jest to drugi po transkrypcji etap realizacji biosyntezy białka (jeden z aspektów przekazywania informacji genetycznej). Translacja jest właściwą syntezą białka i mRNA służy tutaj do przełożenia sekwencji nukleotydów na sekwencje aminokwasów. Aby jednak ten proces mógł zajść musimy zapisać, które aminokwasy kodują odpowiednie nukleotydy.

Zapisujemy nazwy aminokwasów w postaci skrótów 1 i 3 literowych:

calenazwy={"Glycine","LAlanine","LSerine","LProline","LValine","LThreonine","LCysteine","LIsoleucine","LLeucine","LAsparagine","LAsparticAcid","LGlutamine","LLysine","LGlutamicAcid","LMethionine","LHistidine","LPhenylalanine","LArginine","LTyrosine","LTryptophan",stop,"!Unknown"};

skr3={"Gly","Ala","Ser","Pro","Val","Thr","Cys","Ile","Leu","Asn","Asp","Gln","Lys","Glu","Met","His","Phe","Arg","Tyr","Trp", stop,"!??"};

skr1=skr3/.{"Phe"->"F","Tyr"->"Y","Trp"->"W","Arg"->"R","Lys"->"K","Asp"->"D","Glu"->"E","Asn"->"N","Gln"->"Q"};

nazwy1=Join[Table[StringTake[skr1[[i]],1],{i,1,20}],{stop,"?"}];

nazwy=Table[{skr3[[i]],nazwy1[[i]],calenazwy[[i]]},{i,1,22}]

Korzystając z bazy ChemicalData znajdujemy kodon kodujący dany aminokwas:

Table[{amino[[i]],ChemicalData[amino[[i]],"Codons"]},{i,1,20,1}]

kodony={{"GGU","GGC","GGA","GGG"},{"GCU","GCC","GCA","GCG"},{"UCU","UCC","UCA","UCG","AGU","AGC"},{"CCU","CCC","CCA","CCG"},{"GUU","GUC","GUA","GUG"},{"ACU","ACC","ACA","ACG"},{"UGU","UGC"},{"AUU","AUC","AUA"},{"UUA","UUG","CUU","CUC","CUA","CUG"},{"AAU","AAC"},{"GAU","GAC"},{"CAA","CAG"},{"AAA","AAG"},{"GAA","GAG"},{"AUG"},{"CAU","CAC"},{"UUU","UUC"},{"CGU","CGC","CGA","CGG","AGA","AGG"},{"UAU","UAC"},{"UGG"},{"UAA","UAG","UGA"}};

transkrypcja2=Table[First[Flatten[Position[kodony,codons[[i]]]]],{i,1,64}];

Piszemy funkcję, która będzie nam dawać pełną nazwę aminokwasu oraz jego skróty, funkcję tą wykorzystujemy w translacji:

tRNA[rna3_?StringQ,n_]:=Module[{rna3l},
rna3l=StringReplace[rna3,{"CCX"->"CCA","GGX"->"GGA","CGX"->"CGA","GCX"->"GCA","ACX"->"ACA","UCX"->"UCA","CUX"->"CUA","GUX"->"GUA"}];
If[StringFreeQ[rna3l,"X"],nazwy[[transkrypcja2[[First[Flatten[Position[codons,rna3l]]]]],n]],nazwy[[22,n]]]]

Następnie piszemy właściwą funkcję dokonującą translacji:

translacja[transkr_,n_]:=Table[tRNA[transkr[[i]],n],{i,1,Length[transkr]}]

Przykład transkrypcji i translacji

Pobieramy pierwsze (10^5) nukleotydów sekwencji ludzkiego chromosomu 1.

chromosom1=GenomeData[{"Chromosome1",{1,(10^5)}}];

Dokonujemy transkrypcji tej sekwencji przy pomocy napisanych wcześniej funkcji:

transkr=transkrypcja1[chromosom1]

Następnie dokonujemy translacji otrzymanej sekwencji, przy czym ,,Stop" oznacza tzw. ,,kodony stop", które nie kodują żadnych aminokwasów.

transl= translacja[transkr,1]

Ponieważ kodony stop przerywają proces translacji, białko ulega wtedy modyfikacji polegającej na usuwaniu początkowych aminokwasów, a translacja zachodzi na odcinku mRNA zawierajacym jak najmniej kodonów stop.

Replikacja, delecja, mutacja i insercja DNA

Replikacja – jest to duplikacja (powielanie) nici DNA. Jest to proces semikonserwantywny co oznacza, że w każdej uzyskanej nici DNA będzie jedna nić macierzysta i jedna nowa.
Mutacja – każda zmiana w DNA, która nie jest skutkiem procesów rekombinacyjnych.
Delecja – jest to jeden z typów mutacji genowych, polega na utracie jednej lub więcej par nukleotydów z DNA.
Insercja – najczęstsza mutacja genetyczna polegająca na wstawieniu sekwencji DNA do genu.

Piszemy funkcję, która dokonuje replikacji, delecji, mutacji i insercji:

(*Replikacja*)
Replikacja[dna_,pm_,del_,ins_]:=Module[{nukleotydy,randomDNA,PMutacje,len,dna2,k,l,delPos,delMin,delMax,insPos,insMin,insMax,insStrings},
len=StringLength[dna];
(* delecja *)
k=1;l=1;delPos={};delMin=1;delMax=10;
While[k<len+1,
If[RandomChoice[{1-del,del}->{0,1}]==0,k++,
l=Min[k-1+RandomInteger[{delMin,delMax}],len];
AppendTo[delPos,{k,l}];
k=l+1;
]];
dna2=StringReplacePart[dna,"",delPos];
len=StringLength[dna2];

nukleotydy={"A","C","G","T"};
(* mutacje *)
PMutacje[N_]:=If[RandomChoice[{1-pm,pm}->{0,1}]==0,N,RandomChoice[Complement[nukleotydy,{N}]]];
dna2=Table[PMutacje[Characters[dna2][[i]]],{i,len}]//StringJoin;

(* insercje *)
insPos={};insMin=2;insMax=10;insStrings={};
randomDNA[lenR_]:=StringJoin[Table[RandomChoice[nukleotydy],{i,1,lenR}]];
For[k=1,k<(len+1)+1,k++,
If[RandomChoice[{1-ins,ins}->{0,1}]==1,
AppendTo[insPos,{k,k-1}];
AppendTo[insStrings,randomDNA[RandomInteger[{insMin,insMax}]]];
]];
dna2=StringReplacePart[dna2,insStrings,insPos];

dna2
]

Przykład

Nukleotydy= {"A","C","G","T"};

randomDNA[len_]:=StringJoin[Table[RandomChoice[Nukleotydy],{i,1,len}]]

Zakres długości łańcucha DNA

lenRange={50,100};

Losowo ułożona sekwencja DNA

DNA=randomDNA[RandomInteger[lenRange]]

W powyższej sekwencji przeprowadzamy delecje:

DNA2=Replikacja[DNA,.1,.0,.0]

Mutacja:

DNA3=Replikacja[DNA,.0,.05,.0]

Insercja:

DNA4=Replikacja[DNA,.0,.0,.05]

Oraz replikacja:

DNA5=Replikacja[DNA,.1,.05,.05]

Podobieństwa sekwencji

Aby porównać dwie sekwencje DNA możemy zastosować algorytm NeedlemanWunschSimilarity z użyciem macierzy BLOSUM62 (podaje liczbę par nukleotydów) lub algorytm SmithWatermanSimilarity z użyciem macierzy PAM250 (podaje liczbę pokrywajacych się par nukleotydów).

NeedlemanWunschSimilarity

NSWB62[a_?LetterQ,b_?LetterQ]:= NeedlemanWunschSimilarity[a,b,SimilarityRules->"BLOSUM62"]

amino={"A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V"};

len=Length[amino]

amino1= StringJoin[" ",#]&/@amino;

Table[IntegerPart[If[i==j,NSWB62[amino[[i]],amino[[j]]],0]],{i,len},{j,len}];

Table[IntegerPart[If[i>j,NSWB62[amino[[i]],amino[[j]]],0]],{i,len},{j,len}];

%+Transpose[%]+%%//
TableForm[#,TableHeadings->{amino1,amino1},TableAlignments->Right,TableSpacing->{2,1}]&

SmithWatermanSimilarity

SWSP250[a_?LetterQ,b_?LetterQ]:=SmithWatermanSimilarity[a,b,SimilarityRules->"PAM250"]

Table[IntegerPart[If[i==j,SWSP250[amino[[i]],amino[[j]]],0]],{i,len},{j,len}];

Table[IntegerPart[If[i>j,SWSP250[amino[[i]],amino[[j]]],0]],{i,len},{j,len}];

%+Transpose[%]+%%//TableForm[#,TableHeadings->{amino1,amino1},TableAlignments->Right,TableSpacing->{2,1}]&

Optymalne dopasowanie dwóch sekwencji DNA

Polega na takim ustawieniu sekwencji względem siebie, aby miały jak najwiekszą liczbę pokrywających się par nukleotydów.

Podajemy sekwencje, które chcemy ze sobą porównać:

s1="AGCCTTACGAAATATGTAATCGAAGATT"<>"AATATACTACCTGCGCTAGTCGAGCGTC"<>"GAAGAGGAGCAGTGTAAAGCGGAGCGCG"<>"ATGGAATAGATGCGTGGGTGTTCCAAGG"<>"ACAGATGCGGTGTTTGCCCCCCTAC"
s2="ACACACACACACCAACAGCCTTACGTAA"<>"TATGTAATCGAAGATTAATATACCCTAC"

Piszemy funkcję, która będzie przesuwać sekwencje względem siebie pokazując w jakich pozycjach dopasowanie jest największe:

gap[n_]:=Nest[StringJoin["-",#]&,"-",n-1] 
eval[ch1_,ch2_]:=If[ch1==ch2&&ch1≠"-"&&ch2≠"-",1,0]
Eval[s1_,s2_,ali_?NumberQ]:=Module[{S1,S2,len2,s1l,s2l},
If[StringLength[s1]≥StringLength[s2],S1=s1;S2=s2,S1=s2;S2=s1];
len2=StringLength[S2];
s1l=gap[len2]<>S1<>gap[len2]//Characters;
s2l=S2//Characters;
Sum[eval[s1l[[i+ali]],s2l[[i]]],{i,len2}]]

Wykres przedstawia ilość par takich samych nukleotydów w zależności od przesunięcia sekwencji:

EvalT=Table[Eval[s1,s2,k],{k,0,StringLength[s1]+StringLength[s2]}];
ListPlot[%,Joined->True,InterpolationOrder->0,AxesLabel->{"przesunięcie","liczba dopasowań"}]

Best=EvalT//Sort//Last;
pBest=Subtract[#,1]&/@Flatten[Position[EvalT,Best]];

Za pomocą funkcji Manipulate możemy bez problemu zobaczyć w jakim ułożeniu dopasowanie funkcji jest najlepsze:

sh="";
If[StringLength[s1]<StringLength[s2],s1=sh;s1=s2;s2=sh];
len2=StringLength[s2];
s1fix=gap[len2]<>s1<>gap[len2];
s1l=Characters[s1fix];
s2l=Characters[s2];
Manipulate[
{{"Ocena: Suma dopasowań = "<>ToString[Eval[s1,s2,ali]]<>" dla przesunięcia "<>ToString[ali]},{},
{Grid[{
Table[eval[s1l[[i+ali]],s2l[[i]]],{i,len2}],
Table[s1l[[i+ali]],{i,len2}],
Table[s2l[[i]],{i,len2}]
},ItemStyle->{Rule[#,Red]&/@Flatten[Position[Table[eval[s1l[[i+ali]],s2l[[i]]],{i,len2}],1]],1->16}
]
},{},{"Maksymalna suma to "<>ToString[Best]<>" dla przesunięć: "<>ToString[pBest]}
}//TableForm,
{ali,0,StringLength[s1]+len2,1}]

Dopasowanie sekwencji przy użyciu funkcji SequenceAlignment

Od przykładu podanego powyżej ten sposób rożni się tym, iż wstawiamy tzw. przerwy, aby dopasowanie było jak największe. Piszemy funkcje, które pozwolą nam zastosować SequenceAlignment i utworzą przerwy.

Podajemy sekwencje, które chcemy porównać:

s3="TATGGCAGGTTAGACGGAGG";
s4="CTATGGCAGCAGTAGCTGTTAGACGGACG";

{%%,%}//TableForm

eval1[ch3_,ch4_]:=If[ch3==ch4&&ch3≠"-"&&ch4≠"-",1,0]
Eval1[s3_,s4_?StringQ]:=Module[{len,s3l,s4l},
len=Min[StringLength[s3],StringLength[s4]];
s3l=s3//Characters;
s4l=s4//Characters;
Sum[eval1[s3l[[i]],s4l[[i]]],{i,len}]]

gap[n_]:=Nest[StringJoin["-",#]&,"-",n-1]

Funkcja konwertująca format wyjściowy polecenia SequenceAlignment do postaci dwóch sekwencji o takiej samej długości porzez wprowadzenie przerw:

conv[ls_]:=Module[{len1,len2,DeltaLen},
If[StringQ[ls],{ls,ls},
len1=StringLength[ls[[1]]];
len2=StringLength[ls[[2]]];
Which[len1==len2,ls,len2==0,{ls[[1]],gap[len1]},True,{gap[len2],ls[[2]]}]]]

Conv[s3_,s4_,n_]:=Module[{sr},
Which[n==1,sr={{a_,a_}->1, {a_,b_}->-1},n==2,sr={{a_,a_}->0, {a_,b_}->-1},True,sr="PAM250"];
StringJoin/@Transpose[conv/@SequenceAlignment[s3,s4,MergeDifferences->False,SimilarityRules->sr]]]

Conv[s3,s4,1]//TableForm

Funkcja będzie pokazywać również ilość insrecji, delecji, substytucji oraz podobieństwo przy pomocy wspomnianych wsześniej algorytmów:

chngs[s3_,s4_,n_]:=Module[{ls,ident,ins,del},
ls=Conv[s3,s4,n];
{ident=Eval1[ls[[1]],ls[[2]]],ins=StringCount[ls[[1]],"-"],del=StringCount[ls[[2]],"-"],StringLength[ls[[1]]]-ins-del-ident,EditDistance[s3,s4],StringLength[ls[[1]]]-ident,Round[NeedlemanWunschSimilarity[s3,s4]],2*ident-StringLength[ls[[1]]]}]

TableForm[chngs[s3,s4,1],TableHeadings->{{"b.z.","Insercje","Delecje","Substytucje","EditDistance","Ins + Del + Sub","NeedlemanWunschSimilarity","b.z. - (Ins + Del + Sub)"}},TableAlignments->Right]

Podsumowanie

W moim projekcie przedstawiłam przykładowe zastosowanie programu Mathematica do celów bioinformatyki. Ogrom funkcji i możliwości oraz obszerność baz danych pozwala nam dostrzec, iż program ten może mieć zastosowania nie tylko w matematyce i fizyce, ale również w informatycznym podejściu do wybranych problemów biologicznych, choćby jak w podamym przeze mnie bardzo konkretnym przykładzie z dziedziny biochemii.

Bibliografia

Paul G. Higgs, Teresa K. Attwood Bioinformatyka i ewolucja molekularna (Wydawnictwo Naukowe PWN, Warszawa 2008)

http://www.staff.amu.edu.pl/~ewas/pracownia/biomono/

Materiały zdobyte w ramach kursu Bioinformatyka oraz Mathematica I